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SECOND  MOMENT  CLOSURE  MODELING  OF 
COMPLEX  TURBULENT  FLOWS 

AFOSR  GRANT  NUMBER:  FA9550-05-1-0177 
Sharath  S.  Girimaji 
Department  of  Aerospace  Engineering 
Texas  A  &:  M  University,  College  Station,  Texas  77843 

Abstract 

Turbulence  subject  to  unsteady  forcing  can  exhibit  novel  features  that  cannot  be  explained 
using  the  well-known  steady-turbulence  paradigm.  Modeling  and  prediction  of  such  statistically 
unsteady  flows  are  important  in  many  practical  AFOSR  applications:  turbine  flows,  wake-flows 
with  vortex  shedding,  etc.  Further,  many  flow  control  strategies  depend  upon  the  knowledge  of 
unsteady  turbulence  dynamics  to  achieve  the  desired  objectives.  However,  our  understanding 
of  unsteadily-forced  turbulence  dynamics  or  our  ability  to  predict  them  is  inadequate.  Un¬ 
steady  forcing  implies  a  mean  velocity-gradient  field  that  varies  with  time.  Temporally  varying 
mean-velocity  gradient  comprises  of  one  or  more  of  the  following  four  fundamental  forms  of 
unsteadiness:  1)  rotation  of  the  eigen-directions  of  the  mean  strain  rate  tensor  (with  eigen 
values  maintained  constant);  2)  temporal  variation  in  the  eigen  values  of  the  mean  strain  rate 
tensor  (with  eigen  directions  fixed);  3)  rotation  of  the  mean  vorticity-vector  axis  (with  vorticity 
magnitude  maintained  constant);  and  4)  temporal  variation  of  the  mean  vorticity-vector  mag¬ 
nitude  (with  the  direction  fixed).  For  example,  homogeneous  shear  flow  in  a  rotating  reference 
frame,  when  considered  in  an  inertial  frame,  is  a  combination  of  a  mean  strain-rate  tensor 
with  rotating  eigen  vectors  (and  constant  eigen  values)  and  a  constant  vorticity  vector.  In  this 
project  we  undertook  a  systematic  study  of  these  various  forms  of  unsteadiness.  Two  studies 
were  completed  and  the  results  reported  in  two  journal  publications  -  in  Physics  of  Fluids  and 
Journal  of  Fluid  Mechanics.  Two  other  studies  were  initiated.  All  these  efforts  are  detailed  in 
this  final  report. 


1  Introduction 

Over  the  last  several  decades  important  progress  has  been  made  in  our  capability  to  model  arid 
predict  steadily-forced  turbulent  flows  with  nearly  rectilinear  streamlines.  The  next  engineering 
challenge  is  the  modeling  of  turbulent  flows  subject  to  complicating  influences  such  as  unsteady 
forcing  (due  to  external  influences  or  internal  instabilities),  streamline  curvature,  frame-rotation 
and  heat  release.  Many  of  the  applications  of  interest  to  the  Air  Force  involve  flows  in  which 
turbulence  is  strongly  modified  by  the  above  factors.  It  is  well  known  that  such  flows  can  be  best 
modeled  in  mathematically  rigorous  manner  only  at  the  level  of  second  moment  closure.  The 
objective  of  the  present  research  is  two-fold:  (i)  perform  direct  numerical  simulations  (DNS) 
and  rapid  distortion  theory  (RDT)  calculations  of  unsteadily  forced  turbulent  flows  to  infer 
the  behavior  of  pressure-strain  correlation  and  other  key  turbulent  processes;  and  (ii)  develop 
pressure-strain  correlation  and  other  closure  models  which  incorporate  the  requisite  complex 
physics  and  arc  consistent  with  constraints  based  on  the  time-reversibility  characteristics  of  the 
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Navier-Stokes  equations.  These  constraints  will  serve  specific  purposes  for  enhancing  the  model 
capabilities  in  unsteady  flows. 

Over  the  last  two  decades,  direct  numerical  simulations  (DNS)  and  rapid  distortion  theory 
(RDT)  have  played  a  crucial  role  in  providing  important  insight  into  fundamental  turbulence 
processes  leading  to  the  development  of  improved  closure  models.  However,  much  of  this 
progress  is  restricted  to  turbulent  flows  subject  to  steady  forcing.  There  are  numerous  turbulent 
flows  of  interest  to  the  Air  Force  that  are  subject  to  time-dependent,  (unsteady)  forcing:  e.g., 
flows  in  turbines,  internal  combustion  engines,  vortex-shedding  in  wakes.  Furthermore,  there 
are  many  flow  control  strategies  attempting  to  take  advantage  of  unsteady  turbulence  dynamics 
to  modify /control  flow.  Our  understanding  of  fundamental  turbulence  processes  in  flows  with 
rapidly  varying  mean  flows  is  inadequate.  To  accumulate  a  reliable  knowledge  base  on  these 
flows,  we  must  revisit  many  of  the  canonical  flows  and  re-examine  the  turbulence  features  in 
the  context  of  unsteady  forcing. 

Classification  of  unsteadiness.  Detailed  investigations  of  turbulence  subject  to  unsteady 
forcing  arc  relatively  recent.  Unsteady  forcing  implies  a  time- varying  mean- velocity  gradi¬ 
ent  which  encompasses  a  very  wide  range  of  possibilities.  Presumably,  different  types  of  un¬ 
steadiness  can  have  vastly  different  effects  on  turbulence.  This  renders  a  systematic  study  of 
unsteadily-forced  turbulence  rather  difficult.  However,  some  progress  can  be  made  by  recog¬ 
nizing  that  any  arbitrary  unsteadiness  in  the  mean-velocity  gradient  can  be  expressed  as  a 
combination  of  four  elementary  forms  of  variations:  1)  temporal  changes  in  the  eigen- values 
of  the  mean  strain  rate  tensor,  with  eigen-directions  fixed;  2)  rotation  of  the  eigen-directions 
of  the  mean  strain  rate  tensor,  with  eigen-values  maintained  constant;  3)  temporal  changes  of 
the  mean  vorticity-vector  magnitude,  with  the  direction  fixed;  and,  4)  rotation  of  the  mean 
vorticity-vcctor  axis,  with  its  magnitude  maintained  constant.  Our  ultimate  goal  is  to  study, 
in  isolation,  the  effects  of  each  type  of  unsteadiness  on  turbulence. 

Studies  performed  during  funding  period.  An  idealized  flow  which  has  contributed 
greatly  to  our  current  understanding  of  steadily-forced  turbulence  is  the  homogeneous  shear 
flow.  We  revisit  this  flow  and  study  it  in  the  context  of  unsteady  forcing.  Our  objective  is  to 
perform  DNS  and  RDT  of  homogeneous  flows  subject  to  periodic  and  rotating  mean  shear.  The 
goal  is  to  gain  important  insight  into  the  various  physical  processes  in  statistically  unsteady 
turbulent  flows  and  evaluate  the  validity  of  current  turbulent  closures.  We  will  then  attempt  to 
derive  advanced  second-moment  closure  models  that  can  effectively  simulate  these  effects.  It  is 
our  expectation  that  these  second  moment  closure  models  will  serve  as  the  starting  point  from 
which  sophisticated  subgrid  stress  closure  models  can  be  derived  for  high-resolution  computa¬ 
tions.  Two  studies  have  been  completed  and  two  others  initiated  during  this  funding  period. 
We  will  detail  these  studies  later. 

Issues  investigated.  We  address  the  following  issues:  (a)  Effect  of  forcing  frequency  on 
the  evolution  of  turbulence  parameters  such  as  kinetic  energy,  dissipation  and  Reynolds  stress 
anisotropy;  (b)  Budgets  of  Reynolds  stress  and  kinetic  energy  evolution;  (c)  Phase  lag  between 
applied  mean  strain  and  Reynolds  stress;  and  (d)  Comparison  between  DNS  results  and  second 
moment  closure  models. 
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This  project  summary  presents  the  most  prominent  results  reported  in  two  accepted  journal 
articles: 

1.  Yu  and  Girimaji  (2006)  ‘DNS  of  homogeneous  shear  subject  to  periodic  turbulence’.  To 
appear  in  Journal  of  Fluid  Mechanics. 

2.  Girimaji,  O’Neil  and  Yu  (2006)  ‘Rapid  distortion  analysis  of  turbulence  subject  to  rotating 
shear’.  To  appear  in  Physics  of  Fluids. 

The  reader  is  referred  to  these  articles  for  more  detailed  results  and  discussions. 


2  Periodic  shear  flow 

The  first  flow  is  a  homogeneous  periodic  shear  flow  where  the  mean  velocity  gradient  is  given 

by 

d\J ' 

(0  =  SmaxSin (tAj£)<5fi<Jj2-  (1) 

(We  use  the  following  notation:  V,  U  and  u  represent  the  total,  mean  and  fluctuating  velocity 
vectors)  Thus,  the  mean-strain  tensor  has  fixed  eigen  directions  and  temporally  sinusoidal 
eigen  values  (elementary  form  1  in  the  above  list).  In  our  DNS  (direct  numerical  simulation) 
investigation,  the  desired  velocity  gradient  is  produced  by  introducing  a  body  force  which  is 
a  deterministic  function  of  space  and  time.  To  sustain  the  initial  homogeneous  mean  shear,  a 
constant  mean  pressure  gradient  field  is  sufficient. 


Governing  Equations.  To  achieve  time-varying  shear  flow,  additional  body  force,  F (xj), 
is  required.  For  this  purpose  we  will  consider  a  body  force  that  is  a  function  of  space  and 
time,  but  completely  deterministic:  F  does  not  contain  a  fluctuating  part.  Now  the  mean  and 
fluctuating  velocity  equations  can  be  written  as: 


dUt  dUi  d(uiUk) 
dt  kdxk  dxk 


dui  n  dui 
~dt  +  UkdTk  + 


d 

dxk 


(UiUk  -  ( UiUk )) 


dP_ 

dxi 


+  Fi- 


dp'  dUi  d2Ui 

dxi  Uk  dxk  ^  dxkdxk 


(2) 

(3) 


The  fluctuating  velocity  and  pressure  field  equations  are  unaffected  by  the  forcing,  provided 
the  force  is  purely  deterministic.  The  applied  body  force  is 


F  =  UjSmaAy  -  Vref)  COS(( jt)i,  (4) 

where  i  is  the  unit  vector  along  x-dircction.  Also,  ( x.y )  correspond  to  (xi,x2)  directions  and 
yref  is  a  reference  point  which  we  take  to  be  at  the  center  of  the  computational  domain. 
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2.1  Results. 


In  this  section,  we  will  describe  the  observed  results.  Detailed  analysis  and  explanations  arc 
provided  in  the  next  section.  Figure  1  shows  the  evolution  of  the  kinetic  energy  for  various 
forcing  frequencies  over  an  extended  period  of  time.  It  is  clear  that  three  types  of  responses 
are  possible  depending  on  the  frequency  of  forcing.  The  three  regimes  are  identified  as  low-, 
intermediate-  and  high-frequency  regimes. 

When  the  forcing  frequency  is  low  (A/5  <  0.14),  the  evolution  of  the  turbulent  kinetic 
energy  ultimately  asymptotes  to  a  constant  value  which  is  smaller  than  the  initial  value.  At 
intermediate  forcing  frequencies  (A/5  is  between  0.14  and  0.3),  k  increases  at  long  times.  At 
higher  forcing  frequencies(A/5  >  0.3),  the  evolution  of  k  shows  a  periodic  variation,  with  the 
mean  in  the  vicinity  of  initial  value. 

In  Fig.  2,  the  kinetic  energy  evolutions  of  A/5  =  0.16  case  computed  with  three  different 
wavenumber  sets  (of  6078,  13768  and  24318  discrete  wavenumbers)  arc  compared.  The  three 
plots  arc  absolutely  indistinguishable  demonstrating  that  our  results  arc  not  dependent  on 
wave-number  discretization.  All  the  other  results  presented  in  this  study  employ  13768  discrete 
wave-numbers. 

2.2  Low  frequency  Regime 

Figure  3  shows  the  evolution  of  k  in  the  low  frequency  regime.  For  steady  shear  (A/5  =  0), 
k  increases  exponentially.  At  the  lowest  non-zero  frequency  (A/5  =  0.02),  k  increases  initially 
and,  after  reaching  its  maximum  value,  decreases  monotonically  to  an  asymptotic  value.  The 
same  trend  continues  as  frequency  increases,  with  the  initial  peak  and  final  asymptotic  value 
decreasing  progressively.  At  the  highest  frequency  in  this  range,  the  kinetic  energy  is  initially 
periodic  with  cycle-mean  gradually  declining.  At  later  stages,  it  asymptotes  monotonically  to 
a  small  value.  It  appears  that  A/5  =  0  ease  is  a  singular  limit  in  the  low  frequency  regime  as 
the  behavior  at  vanishingly  small  values  is  quite  different. 

Figure  4  shows  the  evolution  of  b\\.  For  steady  shear  (A/5  =  0),  6n  asymptotes  to  the 
maximum  allowable  value  of  2/3.  This  indicates  that  almost  all  of  the  kinetic  energy  in  the 
steady  forcing  case  is  contained  in  the  (u\Ui)  component.  For  non-zero  frequencies,  bn  and  b2 2 
(figure  not  shown)  converge  to  —1/3.  Thus,  all  of  the  kinetic  energy  is  contained  in  (U3U3).  Prior 
to  the  asymptotic  stage,  the  evolution  of  bu  is  different  in  the  A/5  =  0.02  and  0.14  eases.  In 
the  A/5  =  0.02  ease,  61 1  grows  initially  reaching  a  maximum,  and  then  decreases  monotonically 
to  —1/3.  In  the  A/5  =  0.14  ease,  bu  decreases  first,  then  oscillates  with  decreasing  cycle-mean 
and  asymptotes  to  —1/3  at  a  slower  rate. 

In  Figure  5,  for  A/5  =  0  we  see  that,  as  expected,  P\\  is  always  positive  and  grows  with 
time.  This  leads  to  monotonic  growth  of  kinetic  energy  at  long  times.  For  non-zero  frequencies, 
the  evolution  of  Pu  is  quite  different.  The  production  can  be  positive  or  negative  at  early  times, 
but  is  mostly  negative  at  later  times  and  subsequently  asymptotes  to  zero. 

Figure  6  shows  the  evolution  of  the  pressure-strain  redistribution  component  4>u.  In  the 
zero- frequency  (constant  shear)  and  A/5  =  0.02  eases,  the  pressure-strain  correlation  always 
removes  energy  from  (iqui).  In  the  A/5  =  0.02  ease,  there  is  a  period  of  time  (5  <  St  <  10) 
when  both  production  and  pressure-strain  correlation  are  negative.  In  the  A/5  =  0.14  ease, 
$11  always  counters  the  effects  of  Pn.  For  all  frequencies,  asymptotes  to  zero. 
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2.3  Intermediate  frequency  Regime 

The  evolution  of  the  kinetic  energy  for  the  intermediate  range  of  frequencies  (0.14  <  A/5  <  0.3) 
is  shown  in  Fig.  7.  For  all  frequencies,  k  grows  over  long  times.  However,  the  growth  patterns 
are  not  same  for  different  frequencies.  In  the  A/S  =0.16  case,  k  increases  initially  for  a  short 
time  and  decreases,  apparently  asymptoting  to  a  lower  value.  Then,  rather  unexpectedly,  at 
St  =20,  k  increases  dramatically  followed  by  a  decline  that  negates  some,  but  not  all,  of  the 
gain.  Then  it  appears  steady  for  a  period  of  time.  Again,  abruptly,  the  growth-decay-stcady 
pattern  reappears.  In  the  A/S  =0.20  case,  k  varies  periodically,  and  the  cycle-mean  increases 
slowly  (compared  to  the  A/S  =0.16  case)  in  time. 

Figure  8  shows  the  evolution  of  bu.  We  can  sec  that  &n  exhibits  strong  oscillations  at  both 
frequencies.  In  the  A/5  =0.16  case,  the  average  magnitude  of  bn  is  larger  than  that  in  the 
A/S  =0.20  case,  which  reflects  the  different  structure  of  anisotropy  in  the  two  cases.  Both  Pn 
and  <t>u/k  also  exhibit  strong  oscillations  (Figs.  9  and  10).  In  these  cases,  <f>n  always  counters 
the  effects  of  Pu. 

2.4  High  Frequency  Regime 

At  high  frequencies  (A/5  >  0.30),  turbulence  statistics  vary  almost  periodically  as  shown  in 
the  Figs.  11-14.  The  amplitude  of  the  variation  is  smaller  than  those  in  low  and  intermediate 
regime.  Within  this  range,  the  amplitude  decreases  as  the  forcing  frequency  goes  up.  The 
cycle-mean  of  k  is  slightly  larger  than  its  initial  value,  while  the  cycle-mean  of  6U,  Pu  and  4>n 
arc  all  around  zero.  Clearly,  the  turbulence  is  close  to  its  initial  isotropic  state  at  all  times. 
Production  and  pressure-strain  distribution  appear  to  affect  the  velocity  field  very  little. 

2.5  Effects  of  the  phase  of  forcing  shear. 

Thus  far  we  have  discussed  the  effects  of  forcing  (shear)  frequency  on  turbulence.  Now,  we  will 
briefly  discuss  the  effects  of  initial  phase  of  the  periodic  shear.  For  any  arbitrary  initial  phase 
that  a  new  set  of  coordinates  can  be  found  in  which  only  the  (1,2)  component  of  the  mean 
velocity  gradient  tensor  is  non-zero.  Since  the  fluctuating  velocity  field  is  initially  isotropic, 
the  coordinate  transformation  will  not  change  the  statistical  description  of  the  initial  velocity 
fluctuations.  Thus,  all  of  the  preceding  results  pertaining  to  turbulent  kinetic  energy  and  its 
production  will  be  invariant  to  the  initial  phase  of  the  periodic  shear. 

Before  concluding  this  section,  it  is  important  to  point  out  that,  if  the  rapid  pressure 
term  is  omitted  from  the  calculations,  the  turbulent  statistics  vary  purely  periodically  for 
all  frequencies.  Thus,  all  of  the  non-periodic  features  of  the  observed  behavior  at  low  and 
intermediate  frequencies  must  be  attributed  to  the  effects  of  rapid  pressure  term.  In  the  next 
section,  we  will  attempt  to  understand  and  explain  the  role  of  rapid  pressure-strain  correlation. 

2.6  Discussion 

Many  important  aspects  of  rotat.ing-shear  RDT  results  presented  in  the  previous  section  differ 
significantly  from  those  of  homogeneous  shear  flow  in  a  rotating  reference  frame  and  of  periodic 
shear  flow  (mean  velocity  gradient  given  by  Eq.  (1)).  The  most  important  difference  is  at  small 
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unsteadiness  frequencies.  In  the  present  ease,  at  low  forcing  frequencies,  the  kinetic  energy  level 
declines  and  attains  a  steady  state  value  that  is  smaller  than  the  initial  value.  On  the  contrary, 
in  the  other  two  flows,  asymptotic  growth  similar  to  the  steady  forcing  ease  is  observed.  This 
growth  rate,  however,  decreases  with  increasing  unsteadiness  frequency  in  both  homogeneous 
shear  in  rotating  reference  frame  and  periodic  shear  (Eq.  (1))  eases.  Thus,  the  steady-forcing 
behavior  is  a  natural  limiting  case  of  the  vanishing  unsteadiness  frequency  trend  in  those  flows. 
In  the  present  case,  the  steady-forcing  behavior  is  distinctly  different  from  the  small  frequency 
trend.  At  intermediate  unsteadiness  frequencies,  turbulence  grows  in  spurts  in  the  present  ease. 
This  type  of  sporadic  growth  pattern  is  not  seen  in  other  unsteady  flows.  The  physical  reasons 
for  the  growth  and  the  intermittent  manner  in  which  it  happens  arc  unclear.  On  the  other 
hand,  the  behavior  at  high  forcing  frequency  is  quite  understandable.  When  the  shear-rotation 
timescale  is  much  smaller  than  that  of  turbulence  evolution,  production  and  rapid  pressure- 
strain  correlation  process  can  only  affect  small  perturbations  on  the  initial  field.  We  will  now 
attempt  to  provide  a  qualitative  explanation  for  the  observed  apparently  aberrant  behavior  at 
low  and  intermediate  frequencies. 

Steady  shear:  To  explain  the  observed  behavior  at  this  limit,  we  first  recapitulate  the  RDT 
results  of  steady  homogeneous  shear  (see  also  Pope  text  book:  pp.  419-421).  In  the  steady- 
forcing  ease,  all  unit  wavenumber  vectors  evolve  toward  (0,  -1,  0)  indicating  that  all  vectors 
align  along  the  negative  2-dircction  (sec  Eq.  (10)).  Then,  as  a  consequence  of  continuity  u2 
must  vanish.  The  anisotropy  tensor,  at  intermediate  times  (St  <  5),  approaches  the  two- 
componcntial  state  as  (u2u2)  vanishes.  The  (u3u3)  component  does  not  change  significantly  as 
it  experiences  no  production.  On  the  other  hand,  (u\U\)  grows  rapidly.  At  long  times,  turbulent 
fluctuations  approach  a  one-componential  state  as  the  magnitude  of  (iqui)  vastly  exceeds  that 
of  (u3u3).  Thus,  in  general,  the  dominant  effect  of  steady  shear  is  to  augment  Reynolds  stress 
component  in  the  mean-flow  direction  ((uiiq)  in  this  ease)  and  suppress  the  component  in  the 
normal  direction  ((ii2U2)).  The  augmentation  of  the  streamwise  component  is  due  to  the  effect 
of  production  and  suppression  of  the  flow-normal  component  is  due  to  the  action  of  pressure- 
strain  correlation.  We  will  now  argue  that  this  suppression  of  the  flow-normal  component  is 
key  to  explaining  the  observed  behavior. 

Low  forcing  frequencies:  Now  let  us  consider  the  slowly  rotating  shear  case.  Initially,  the 
orientation  of  shear  is  such  that  (iqi^)  grows  as  seen  from  Fig.  4  and  (u2u2)  rapidly  decreases 
(figure  not  shown).  The  decrease  in  (u2u2)  is  due  to  two  factors.  The  negative  production  (that 
develops  as  shear  rotates)  contributes  to  some  of  the  reduction  and  the  rest  is  due  to  rapid 
pressure  term  as  in  the  ease  of  steady  shear.  During  the  intermediate  stages  of  evolution,  the 
shear  rotates  (at  a  rate  dictated  by  A)  to  an  orientation  so  that  the  velocity  gradient  tensor  is 
of  the  form 


dUi 

dxj 


=  6l2S 


(5) 


The  new  orientation  favors  the  growth  of  (u2u2)  and  reduction  of  (uiiq).  If  (u2u2)  is  driven  to 
zero  in  the  early  stages  before  favorable  shear  alignment  can  augment  it,  then  this  component 
will  remain  zero  at  all  subsequent  times.  This  is  due  to  the  fact  that  favorable  shear  requires  non¬ 
zero  u2  perturbation  to  produce  (u2u2)  component  of  kinetic  energy.  During  this  intermediate 
stage,  the  (iqui)  component  which  had  previously  experienced  growth  is  quickly  depleted  by 
the  action  of  pressure  strain  correlation.  In  the  lowest  frequency  case,  it  is  indeed  driven  to 
zero  during  this  stage  as  seen  from  Fig.  4.  At  slightly  higher  frequencies,  the  depletion  is 
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not  complete  in  one  cycle  of  shear  rotation.  If  the  alignment  of  shear  changes  before  (u\U\)  is 
completely  depleted  in  the  second  quarter  of  the  cycle,  this  component  of  stress  will  increase 
as  it  experiences  favorable  shear  in  the  third  quarter  of  rotation.  Then,  ( U\U\ )  will  alternately 
rise  and  fall,  with  depletion  exceeding  growth  in  each  quarter  cycle.  Ultimately,  the  stress  will 
vanish  over  several  cycles  of  shear  rotation  as  seen  in  A/5  =  0.14  case  in  Fig.  4. 

The  role  of  rapid-pressure  strain  correlation  can  now  be  summarized  as  follows.  This  term 
depletes  the  Reynolds  stress  component,  in  the  flow  normal  direction:  (u2u2)  during  the  first, 
and  third  quarter  and  ( U\U\ )  during  the  second  and  fourth  quarter  of  each  shear  rotation  cycle. 
At  low  frequencies,  the  depletion  is  more  complete  than  at  high  frequencies.  It  is  this  effect 
that  leads  to  the  observed  behavior.  Thus,  in  the  absence  of  pressure-strain  correlation,  the 
Reynolds  stress  behavior  will  be  quite  different. 

Intermediate  forcing  frequencies:  As  the  frequency  of  rotation  increases  further  (0.14  < 
A/5  <  0.16),  more  of  the  stresses  survive  the  depletion  cycle  due  to  rapid  changes  in  the 
orientation  of  strain  and  consequently  the  action  of  pressure.  Once  in  every  few  cycles,  the 
shear  aligns  favorably  with  the  existing  fluctuations  leading  to  the  spurts  of  rapid  growth. 
Following  these  periods  of  fast  growth,  the  changing  strain  orientation  leads  to  fairly  steep 
decline  due  to  the  combined  action  of  pressure-strain  correlation  and  negative  production.  But 
the  decline  is  not  large  enough  to  offset  the  growth,  resulting  in  a  net  gain  of  kinetic  energy  after 
each  growth-decline  cycle.  With  increasing  frequency,  the  growth  spurts  occur  more  frequently. 
However,  the  magnitude  of  growth  diminishes  as  the  requisite  favorable  alignment  occurs  for  a 
shorter  duration. 

High  forcing  frequencies:  When  the  forcing  frequency  is  very  high,  the  strain  rate  rotates 
too  rapidly  and  production  is  not  sustained  for  a  long  enough  period  for  ( U\U\ )  to  grow.  The 
rapid  rotation  also  causes  the  pressure-strain  correlation  to  change  too  quickly  to  suppress 
(u2?x2)  to  any  significant  degree.  As  a  result,  the  kinetic  energy  changes  very  little  from  its 
initial  value. 

Phase  difference  between  mean  strain  tensor  and  Reynolds  stress.  The  evolution 
equation  for  kinetic  energy  is 

dk 

—  =  —2  kbijSij.  (6) 

The  orientation  between  and  Stj  is  crucial  in  determining  the  evolution  of  k.  It  is  the 
action  of  rapid  pressure  that  creates  the  phase  lag  between  the  applied  strain  and  Reynolds 
stress.  Clearly,  the  phase  lag  will  depend  strongly  on  the  shear  rotation  rate.  For  a  given  A/5, 
calculations  show  that  the  phase  differences  between  &n  and  S\\  vary  about  a  median  value 
(a)  as  a  function  of  time.  It  is  this  median  value  that  is  important  for  explaining  the  observed 
behavior  as  a  function  of  forcing  frequency.  It  is  hard  to  directly  compute  this  value  from  the 
data  and  we  will  surmise  this  indirectly  from  the  sign  and  magnitude  of  production.  Positive 
production  will  indicate  favorable  alignment  (n/2  <  a  <  n).  Negative  production  will  indicate 
unfavorable  alignment  (n/2  >  a  >  0). 

1.  In  the  constant  shear  ease,  the  monotonic  growth  in  production  (Fig.  5)  clearly  indicates 
that  the  alignment  between  strain  and  stress  is  favorable  (o'  ~  n). 

2.  In  the  low  non-zero  frequency  eases,  the  production  is  negative  at  long  times  indicating 
that  the  alignment  is  unfavorable  (a  <  7r/2) ,  consistent  with  the  explanation  provided 
earlier  in  this  section. 
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3.  The  nearly  periodic  behavior  of  production  at  the  intermediate  frequencies  (Fig.  9)  shows 
that  the  phase  difference  is  nearly  7t/2  but  with  occasional  favorable  alignment  (a  >  7t/2) 
resulting  in  positive  spikes  in  production  as  explained  earlier  in  the  section.  While  some 
negative  excursions  of  production  arc  also  seen,  the  positive  spurts  outnumber  the  negative 
ones  leading  to  overall  positive  production  during  these  periods.  On  an  average,  in  the 
intermediate  regime,  the  alignment  can  be  taken  to  be  slightly  favorable  (a  >  7r/2). 

4.  At  the  highest  frequencies,  it  is  clear  that  the  phase  difference  is  a  %  7t/2  on  an  average, 
as  no  particular  alignment  is  sustained  for  long  enough  time. 

Simple  mathematical  model:  Now  we  will  attempt  to  construct  a  simple  mathematical 
model  that  can  mimic  the  observed  range  of  behavior.  Consider  the  scalar  equation 

dk 

—  =  —  k  cos(t)cos(t  +  a)  (7) 

In  the  equation,  a  is  the  phase  lag  and  cos(2)  and  cos(t  H-  a)  model  bn  and  Sn,  respectively. 
Integrating  the  above  model,  we  have  the  following  evolution  equation  for  the  kinetic  energy: 

—  e~\t  cos(a)-±sin(2t+a)+±sin(ar)  (g) 

Figure  15  shows  the  evolution  of  k  using  Eq.  (8)  for  several  values  of  a.  We  can  clearly  see 
that  there  are  three  different  regimes  of  responses  depending  on  value  of  a.  For  n/2  <  a  <  n, 
modcl-fc  evolution  is  very  similar  to  that  seen  in  the  intermediate  forcing  regime.  For  a=  7r/2, 
k  is  periodic  and  the  RDT  behavior  in  high  frequency  regime  is  recovered.  When  a  <  n/2, 
k  decays  to  zero.  In  the  RDT  ease,  at  low  forcing  frequencies,  k  settles  down  to  a  constant 
value  lower  than  unity.  Thus,  the  turbulence  responses  to  rotating  shear  can  be  explained 
qualitatively. 

Summary  of  Study  1.  The  main  findings  are: 

1.  At  low  forcing  frequencies,  kinetic  energy  grows.  At  high  frequencies,  k  decays.  The 
critical  frequency  at  which  the  asymptotic  behavior  changes  from  growth  to  decay  is  found 
to  be  around  uj/Smax= 0.5.  At  very  high  frequency  (u;/5max  =  10),  the  periodically-forced 
turbulence  behaves  similar  to  decaying  isotropic  turbulence. 

2.  The  observed  behavior  can  be  explained  in  terms  of  the  stress-strain  phase-lag  <p  at  dif¬ 
ferent  forcing  frequencies.  The  phase-lag  goes  from  n  at  low  frequencies  to  n/2  at  high 
frequencies.  Phase-lag  of  n/2  corresponds  to  no  net  production.  Thus  the  production  de¬ 
creases  progressively  with  increasing  frequency.  When  the  level  of  production  falls  below 
dissipation,  turbulence  cannot  be  sustained. 

3.  Normal  anisotropy  oscillates  about  non-zero  asymptotic  values  in  growth  eases.  In  decay 
eases,  the  cycle-average  of  anisotropy  appears  to  vanish  at  long  time. 

4.  RDT  docs  not  capture  the  frequency-dependence  of  the  asymptotic  behavior.  It  predicts 
periodic  A;-bcha.vior  for  all  frequencies. 
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5.  Second  moment  closures  do  capture  the  asymptotic  frequency  dependence,  although  the 
predicted  critical  value  is  smaller  than  that  observed  in  DNS. 

6.  Comparison  of  RDT  and  DNS  results  appears  to  indicate  that  the  onus  is  on  the  closure 
models  for  the  non-linear  terms  to  produce  the  observed  asymptotic  frequency  depen¬ 
dence. 

The  next  step  is  to  develop  second  moment  closures  that  do  capture  the  behavior  accurately. 
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Figure  1:  Extended  kinetic  energy  evolution 


Figure  2:  Kinetic  energy  evolution  at  various  wavenumber  discretizations:  A/S  =  0.16 
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Figure  3:  Low-frequency  evolution  of  kinetic  energy 
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Figure  4:  Low-frequency  evolution  of  6n 
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Figure  5:  Low-frcqucncy  evolution  of  Pn 
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Figure  6:  Low-frequency  evolution  of  3>n 


Figure  7:  Intermediate-frequency  evolution  of  kinetic  energy 
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Figure  8:  Inter  mediate-frequency  evolution  of  6n 
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Figure  9:  Intermediate-frequency  evolution  of  Pn 
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Figure  10:  Intermediate-frequency  evolution  of  ^E>n 
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Figure  11:  High-frequency  evolution  of  kinetic  energy 
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Figure  12:  High-frequency  evolution  of  b\\ 
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Figure  13:  High-frequency  evolution  of  P\\ 
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Figure  14:  High-frequency  evolution  of  4>n 
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Figure  15:  Evolution  of  k  for  various  values  of  phase  lag,  a 
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3  Rotating  Shear  Flow 


The  second  form  of  unsteady  homogeneous  shear  which  is  being  investigated  in  a  concurrent 
study  is  a  rotating  shear  flow: 


=  5 


sin(2nAt)cos(2nAt)  cos2(2nAt)  0 

—sin2(2nAt)  —sin(27TAt)cos(2nAt)  0 
0  0  0 


(9) 


In  this  flow,  the  eigen  values  of  the  mean-strain  rate  tensor  are  constant  but  the  eigen-directions 
rotate  at  a  constant  rate  (elementary  form  2).  When  the  rotation  rate  is  zero,  the  forcing  reduces 
to  a  simple  homogeneous  shear.  Again,  the  required  time- varying  velocity  field  is  produced  using 
a  body  force.  Due  to  the  basic  difficulties  in  performing  DNS  of  this  flow,  we  will  resort  to  RDT 
computations  and  restrict  our  considerations  to  linear  processes.  This  type  of  unsteadiness  is 
interesting  for  two  reasons:  (i)  it  is  one  of  the  four  elementary  forms  of  unsteadiness  as  discussed 
in  the  introduction;  and  (ii)  the  observed  behavior  as  a  function  of  frequency  is  remarkably 
different  from  that  of  other  unsteady  flows  (periodic  shear  flows  and  homogeneous  shear  flow 
in  rotating  reference  frame).  Rotating  homogeneous  shear  flow  shares  some  common  features 
with  the  homogeneous  turbulence  subject  to  periodic  shear  (elementary  form  1)  and  rotating 
shear  (elementary  form  2).  When  considered  in  the  inertial  frame,  it  is  clear  that  the  eigen 
directions  of  the  applied  shear  rotate  and  turbulence  is  subjected  to  time-dependent  forcing. 
Rotating  homogeneous  shear  flow  has  been  the  subject  of  many  theoretical  and  computational 
studies.  The  major  findings  from  these  studies  are:  (a)  for  small  rates  of  rotation  (D)  relative 
to  shear  S,  turbulence  grows  in  time,  although  slower  than  in  the  unrotated  shear  case;  and, 
(b)  for  high  rates  of  rotation,  turbulent  fluctuations  become  two-dimensional  (as  per  Taylor- 
Proudman  theorem)  and  turbulence  decays  rapidly.  It  is  found  that  the  observed  behavior  is 
well  predicted  by  linear  stability  theory.  Linear  stability  equations  predict  that  turbulence  can 
be  sustained  only  for  the  range  of  parameter  values  0  <  |fi|/S  <  0.5.  This  range  is  very  close  to 
that  observed  in  experiments  and  second-moment  closure  model  calculations.  We  will  examine 
if  turbulence  subject  to  rotating  shear  in  an  inertial  reference  frame  exhibits  similar  behavior. 


Governing  Equations.  The  governing  equations  arc  most  expeditiously  solved  in  Fourier 
space  (/c).  The  linearity  of  the  governing  equation  at  the  RD  limit  permits  the  analysis  of 
velocity  (u(£,t))  and  pressure  (p(k,t))  of  each  Fourier  mode  independently.  The  analysis  is 
typically  performed  in  a  coordinate  frame  moving  with  the  mean  velocity.  Thus  the  wave- 
number  vector  of  a  Fourier  mode  itself  is  a  function  of  time,  Pope  (2000).  It  can  be  shown  that 
the  evolution  equations  for  k( t)  and  u(t)  are  given  by: 


dki  _  ^  dUj 

dt  J  dxi 


(10) 


duj 

dt 


(11) 


Note  that  these  Fourier  coefficients  are  further  subject  to  the  following  incompressibility  con¬ 
dition: 


u j/Sj  =  0. 


(12) 
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The  behavior  of  the  velocity  and  pressure  fields  can  be  determined  by  simple  integration  over 
all  wavenumbers.  Summing  this  equation  over  all  wavenumbers  gives  the  evolution  equation 
for  the  Reynolds  stress  in  the  RD  limit: 

—  (utUj)  =  -(ujUk)  ^  -  (uiUk)  ^  +  Mikjl)  (13) 

where  M^i  is  a  fourth  order  tensor  defined  by 

Mijki  =  ] 'T  UiUjikei.  (14) 

Turbulence  production,  is  identified  as 

,  ,  dUi  ,  dUj  ,  . 

Pii  =  -<»,«*>  Wk  -  <W  05) 

and  the  rapid  pressure-rate-of-strain  tensor  is  defined  as 

=  2^(Mki«  +  Mikjl).  (16) 

Then  Eq.  (13)  can  be  rewritten  as 

|  (ulUj)  =  ^  +  <t>y.  (17) 

Thus,  at  the  rapid  distortion  limit,  Reynolds  stress  evolution  is  governed  entirely  by  the  linear 
processes  of  production  and  rapid  pressure-strain  correlation.  The  first  step  in  performing  the 
RDT  simulation  is  to  initialize  the  turbulence  velocity  field.  In  the  present  study,  we  consider 
an  initially  isotropic  velocity  field  that  is  subsequently  subjected  to  unsteady  forcing. 

3.1  DNS  Results. 

The  effects  of  six  different  forcing  frequencies  on  turbulence  evolution  are  studied:  uj/Smax 
=  0.125,  0.25,  0.50,  0.75,  1.0  and  10.  In  all  the  cases,  the  initial  turbulence  field  is  isotropic 
and  prepared  as  described  earlier  and  allowed  to  evolve.  The  various  statistical  moments 
arc  computed  as  a  function  of  time.  The  results  are  sorted  into  several  different  categories: 
kinetic  energy  and  dissipation;  production;  shear  anisotropy;  normal  anisotropy;  Reynolds 
stress  budget;  and,  limiting  high-frequency  behavior.  Each  category  is  examined  individually 
and  in  relation  with  other  categories. 

3.1.1  Evolution  of  k  and  e 

In  all  the  cases  studied,  turbulence  decays  initially  as  the  production  is  nearly  zero  (due  to 
Reynolds  stress  isotropy)  and  dissipation  is  comparatively  large.  After  the  initial  decay,  subse¬ 
quent  evolution  of  k  is  different  for  different  frequencies  as  shown  in  Fig.  24.  In  u;/5max= 0.125 
and  0.25  cases,  k  grows  rapidly  after  the  initial  period  of  decay.  The  growth,  however,  is  not 
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monotonic  and  cyclic  variations  from  the  mean  trend  is  very  much  in  evidence.  In  the  ease 
of  cj/5=0.5.  the  initial  decay  period  is  followed  by  a  long  duration  of  purely  periodic  behav¬ 
ior.  Beyond  Smaxt  ~  25,  the  kinetic  energy  shows  signs  of  very  slow  growth.  In  u;/5max=0.75 
and  1.0  eases,  k  appears  to  decay  with  each  cycle  of  applied  strain.  It  is  clear  that  two  dis¬ 
tinct  turbulence  responses  arc  possible  depending  on  the  frequency  of  forcing.  At  low  shear 
frequencies,  turbulence  grows  at  long  times,  although  not  necessarily  monotonically.  At  high 
shear  frequencies,  turbulence  decays.  The  switch  from  turbulence  decay  to  growth  occurs  at  a 
frequency  of  approximately  u;/5max  =  0.5.  It  is  likely  that  at  some  normalized  frequency  close 
to  0.5  turbulence  neither  decays  nor  grows.  It  is  however,  difficult  to  precisely  pin-point  this 
value  of  frequency  for  which  the  turbulence  level  remains  close  to  its  initial  condition.  While 
the  turbulence  dynamics  in  the  growth  eases  is  likely  to  be  similar  to  the  steady  forcing  ease, 
the  physics  of  the  decaying  eases  need  further  investigation. 

Two  decaying  eases  (c d/Smax  =  0.5  and  1.0)  arc  studied  in  more  detail  next  in  Figs.  25 
and  26.  In  these  figures,  the  cycles  of  k  variation  are  compared  to  those  of  5-variation  and 
the  behavior  of  dissipation  is  also  examined.  Clearly,  the  frequency  of  variation  k  is  twice 
that  of  shear.  We  identify  five  points  in  time  (A,  B,  C,  D,  and  E)  in  one  period  (T0)  of  shear 
cycle.  Points  A,  C,  and  E  correspond  to  consecutive  zero-crossing  times  of  5  and  points  B  and 
D  correspond  the  consecutive  minima  times  in  k  variation.  In  u;/5max= 1.0  case,  consecutive 
k  peaks  have  progressively  smaller  magnitudes  indicating  decay.  In  the  u/Smax  =  0.5  ease, 
consecutive  peaks  in  the  intermediate  period  of  evolution  arc  nearly  identical. 

In  all  the  eases  considered,  the  evolution  cycles  of  k  and  s  arc  in  phase.  This  observation 
has  important  modeling  implications.  The  current  closure  strategy  of  modeling  production  of 
dissipations  in  terms  of  production  of  kinetic  energy  appears  to  be  valid  for  time-dependent 
forcing  as  well,  at  least  at  these  low  Reynolds  numbers. 


3.1.2  Time  dependence  of  P/s 

A  simple  examination  of  the  kinetic  energy  equation  in  homogeneous  flows 


e 


(18) 


reveals  that  /c-evolution  is  completely  dictated  by  the  difference  between  production  and  dissi¬ 
pation.  In  order  to  understand  the  observed  evolution  of  k,  we  investigate  the  behavior  of  P/s 
ratio.  This  ratio  must  exceed  unity  for  turbulence  to  grow. 

Figures  27  and  28  show  the  evolution  of  P/s.  The  evolution  is  generally  oscillatory  with 
periods  of  negative  production.  Negative  production  is  rarely  seen  in  turbulence  subject  to  con¬ 
stant  forcing  but  may  be  important  in  many  practical  flows.  Over  the  first  half-cycle  of  forcing, 
strain  and  production  arc  in  phase  and  their  frequencies  appear  to  be  the  same.  However,  a 
phase  difference  soon  develops  and  the  frequency  of  the  ratio  ( P/s )  quickly  becomes  twice  that 
of  the  applied  strain  as  can  be  seen  from  monitoring  the  five  reference  points  A,  B,  C,  D  and 
E  (defined  in  Figs.  25  and  26).  At  all  these  points,  production  goes  to  zero.  The  zeros  at  A, 
C  and  E  are  due  to  the  strain-rate  going  to  zero.  The  zeros  at  B  and  D,  it  will  be  seen  soon, 
arc  due  to  vanishing  Reynolds  stress  component  (uiu2).  The  frequency  of  the  production  cycle 
drives  the  kinetic  energy  cycle.  Thus,  the  kinetic  energy  cycle  is  also  twice  that  of  the  applied 
strain  rate  as  seen  previously.  The  cycle-averaged  ratio  grows  gradually  beyond  unity  in  the  low 
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frequency  eases:  c o/Smnx  <  0.5.  At  higher  frequencies,  the  cycle-averaged  ratio  lingers  around 
zero  (c.g.,  uj/Smax  =  1  case  shown  in  Fig.  28).  Thus,  the  observed  k  behavior  can  be  completely 
explained  in  terms  of  the  asymptotic  trends  of  mean  production-to-dissipation  ratio. 

3.1.3  Evolution  of  612 

To  gain  further  insight  into  the  behavior  of  P/e,  we  next  examine  the  evolution  of  shear 
anisotropy  b12,  which  in  conjunction  with  S(t)  determines  production.  In  constant  shear  tur¬ 
bulence,  b\2  and  S  arc  of  opposite  signs. 

The  evolution  of  anisotropy  component  612  in  the  uj/Smax= 0.5  case  is  shown  in  Fig.  29. 
The  time  variation  of  shear  is  also  given  in  the  figure.  For  the  first  quarter  of  the  cycle,  612 
is  in  phase  with  S  and  evolution  is  as  expected.  As  the  value  of  the  applied  strain  begins  to 
diminish  after  peaking,  a  phase  lag  develops.  When  the  strain  goes  to  zero  at  the  half-cycle 
mark,  the  stress  is  clearly  non-zero  highlighting  the  hysteresis  effect.  For  a  substantial  portion 
of  the  third  quarter  of  the  first  cycle,  stress  and  strain  are  of  the  same  sign,  implying  that  the 
production  is  now  negative.  Toward  the  end  of  the  third  quarter,  the  stress  passes  through 
zero  and  production  becomes  positive  again.  Throughout  the  fourth  quarter  of  the  first  cycle, 
the  stress  and  strain  are  of  opposite  signs.  At  the  end  of  the  first  cycle,  the  stress  does  not 
return  to  zero.  The  phase-lag  between  stress  and  strain  continues  to  grow  for  another  cycle 
of  applied  strain.  By  the  end  of  the  second  cycle,  stress  and  strain  lock  into  a  constant  phase 
difference  and  evolve  at  the  same  frequency.  In  this  asymptotic  state,  the  dynamics  is  again 
examined  using  the  reference  points  A,  B,  C,  D  and  E.  We  divide  the  A-E  time  period  into 
four  intervals:  AB,  BC,  CD,  and  DE,  as  shown  in  Fig.  29.  Just  before  time  A,  the  shear  is 
negative,  bn  is  positive,  and  production  is  positive.  In  the  period  AB,  shear  becomes  positive 
while  bn  still  maintains  positive  values.  In  this  period,  production  is  negative  as  shown  in  Fig. 
27.  The  negative  production  is  also  reflected  in  the  rapid  fall  of  k  during  the  same  period  in 
Fig.  25.  It  can  be  seen  in  Fig.  29  that  shear  and  6J2  are  relatively  small  in  period  AB  hence 
the  negative  production  is  relatively  small  in  magnitude.  After  time  B,  production  returns 
to  positive  values  until  time  C.  As  seen  in  Fig.  29,  concurrence  of  large  shear  and  bn  values 
produces  large  positive  production  during  period  BC.  During  period  CD  and  DE,  we  sec  the 
same  kind  of  production  behavior  as  in  periods  AB  and  BC. 

Figure  30  shows  the  evolution  of  bn  for  cj/Smai=1.0  case.  The  differences  between  the 
two  cases  are  clearly  evident.  In  u/Smax  =  1.0  case,  b\2  is  not  initially  symmetric  about  zero 
with  the  negative  values  being  much  larger  than  positive  values.  This  is  caused  by  the  initial 
sign  of  the  applied  shear.  In  our  simulation,  the  initial  shear  is  positive  leading  to  large  initial 
negative  values  of  6]2  in  the  first  half  of  the  first  cycle.  At  the  second  half  of  first  cycle,  the 
shear  is  negative  and  612  is  changing  to  positive  values  to  keep  in  pace  with  the  negative  shear. 
However,  before  the  6i2  can  reach  its  full  potential  magnitude  on  the  positive  side,  the  first  cycle 
rapidly  comes  to  an  end  and  the  applied  shear  is  positive  again.  Thus,  the  positive  excursions 
of  b\2  are  not  as  large  initially  as  the  negative  excursions.  With  the  passage  of  time,  the  bias 
toward  negative  values  decrease  as  the  effect  of  the  initial  sign  of  applied  shear  fades  and  the 
oscillations  become  more  symmetric.  Calculations  with  initial  negative  shear  show  the  opposite 
initial  bias  but  do  not  change  the  observed  asymptotic  behavior  of  the  evolution  of  kinetic 
energy,  production  and  dissipation. 

The  difference  between  the  low  and  high  frequency  cases  at  long  times  is  more  important. 
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Examining  the  time  period  A-E  in  the  two  eases,  it  can  be  seen  that  in  cj/Smax=1.0  ease,  period 
CD  is  much  longer  and  period  DE  is  much  shorter  than  they  arc  in  u;/5max=0.5  ease.  This 
generates  a  large  negative  production  during  time  CD  and  small  positive  production  during 
time  DE  in  uj/Smax=  1.0  ease.  Overall,  the  net  production  during  a  cycle  can  be  judged  by  the 
lengths  of  AB+CD  and  BC+DE  which  correspond  to  the  times  when  production  has  negative 
and  positive  values,  respectively.  If  BC+DE  is  larger  than  AB+CD,  net  production  will  be 
positive.  Otherwise,  the  cycle  average  of  production  will  be  negative.  As  b\2  and  S  vary 
periodically,  these  lengths  can  be  characterized  by  a  single  parameter:  asymptotic  phase  lag 
((f))  between  S  and  fr12.  If  7r/2  <  (f)  <  n,  we  will  have  (BC+DE)>  (AB+CD),  and  the  net 
production  will  be  positive.  If,  on  the  other  hand,  0  <  0  <  7r/2,  we  will  have  (BC+DE)  < 
(AB+CD),  and  the  net  production  will  be  negative.  If  (f>  =  7r/2,  then  (BC+DE)  =  (AB+CD), 
resulting  in  no  net  production  in  a  cycle  of  applied  strain.  For  constant  shear  case,  considering 
a  constant  shear  as  a  pulse  wave  with  infinite  period,  we  have  0=7r.  Thus  the  production  will 
be  always  positive.  In  our  computations,  the  phase  difference  between  stress  and  strain  vary 
initially,  but  lock  into  a  constant  value  at  latter  times.  The  dependence  of  the  lock-in  phase-lag 
as  a  function  of  the  applied  frequency  is  shown  in  the  original  paper.  It  is  seen  that  phase 
lag  goes  from  nearly  n  in  the  very  low  frequency  case  to  about  n/2  in  high  frequency  eases. 
Thus,  the  net  production  decreases  with  increasing  frequency  of  forcing.  It  is  curious  to  note 
that  (p  docs  not  go  below  7r/2  implying  that  net  production  is  always  non-negative.  Thus,  the 
net  energy  transfer  in  each  cycle  is  from  the  mean  to  fluctuating  field,  except  at  very  high 
frequencies  at  which  the  transfer  goes  to  zero. 

In  order  to  sustain  turbulence,  the  production  must  not  only  be  positive,  but  must  also 
exceed  dissipation.  From  our  computations,  it  appears  that  for  u/Smax  >  0.5,  net  cycle  pro¬ 
duction  is  smaller  than  net  cycle  dissipation.  We  identify  the  critical  frequency  as  u;cr  =  0.5^^, 
beyond  which  turbulence  cannot  be  sustained.  At  higher  frequencies,  dissipation  rate  exceeds 
the  rate  at  which  energy  is  transfered  from  mean  to  fluctuating  field. 

3.1.4  Normal  anisotropy 

We  now  turn  our  attention  to  the  normal  components  of  Reynolds  stress  tensor.  In  inhomo¬ 
geneous  flows,  these  components  play  an  important  role  in  determining  the  secondary  flow 
structures  of  the  mean  velocity  field.  The  details  of  the  evolution  of  the  diagonal  anisotropy 
components  in  u;/Smax=0.5,  0.65  and  1.0  cases  are  given  in  Figs.  31,  32  and  33.  The  major 
observations  and  their  implications  are  now  summarized. 

The  frequency  of  all  normal  (or  diagonal)  anisotropies  is  the  same  as  that  of  kinetic  energy 
and  twice  that  of  the  applied  shear  or  off-diagonal  (&12)  anisotropy.  This  is  to  be  expected  as 
the  frequency  of  kinetic  energy  is  also  twice  that  of  shear  due  to  the  influence  of  production. 
That  the  diagonal  Reynolds  stress  component  should  have  twice  the  frequency  off-diagonal 
component  can  also  be  understood  from  a  second  point  of  view.  Considering  one  cycle  of  the 
evolution  of  (u\U2),  we  can  expect  ( UiU\ )  and  (u2u2),  which  arc  always  positive,  to  attain  their 
maximum  values  during  the  maximum  positive  or  maximum  negative  value  of  (u\u2).  The 
diagonal  Reynolds  stress  will  attain  their  minimum  value  when  (u\u2)  crosses  zero  on  its  way 
from  maximum  to  minimum  value.  Thus,  in  the  time  that  it  takes  the  off-diagonal  component 
to  complete  a  half-cycle  from  maximum  to  minimum  through  zero,  the  diagonal  components 
complete  a  full  cycle.  The  difference  in  frequency  between  the  applied  strain  and  diagonal 


23 


stresses  also  has  important  modeling  implications.  Approximations  which  imply  that  stress 
and  strain  have  the  same  frequency  will  be  invalid. 

In  the  kinetic  energy  growth  eases,  the  asymptotic  cycle-mean  anisotropy  values  arc  non¬ 
zero.  For  example,  in  the  case  of  u/Smax  =  0.5,  the  asymptotic  cycle-mean  values  of  bu  and 
622  arc  about  0.15  and  -0.1  respectively.  In  general,  at  all  forcing  frequencies,  component 
6n  oscillates  about  a  positive  mean  value  and  622  about  a  negative  value.  The  magnitude  of 
the  cycle-mean  values  of  diagonal  anisotropies  decrease  with  increasing  frequency.  Thus,  with 
increasing  frequency,  the  diagonal  stresses  tend  closer  toward  isotropic  state.  In  all  the  decay 
eases,  normal  anisotropy  tends  to  zero  asymptotically,  even  if  it  develops  non-zero  values  during 
first  cycle  of  shear.  Thus,  the  decay  cases  evolve  asymptotically  to  isotropic  turbulence. 

3.1.5  Budget  of  Reynolds  stress 

Reynolds  stress  budgets  offer  more  detailed  insights  into  the  physics  of  turbulence  processes  and 
arc  very  important  for  high-order  moment  closure  model  development.  We  will  now  investigate 
the  budgets  of  the  various  Reynolds  stress  components  in  the  u;/5mai  =  1.0  case.  Although  the 
computations  are  performed  with  the  lattice  Boltzmann  equation,  the  budgets  will  be  discussed 
in  the  context  of  the  Navicr-Stokcs  equations.  (The  Reynolds  stress  budget  equations  are  given 
in  Section  3.)  Following  Pope  (2000),  production,  pressure-strain  correlation  and  dissipation 
are  represented  by  P^,  Rij  and  £{j.  Wherever  possible,  these  budgets  will  be  put  in  perspective 
with  the  budgets  in  constant  shear  homogeneous  turbulence  (Rogers  1986). 

First  and  foremost,  we  must  ensure  that  the  numerical  data  is  of  high  enough  numerical 
fidelity  that  the  various  high-order  derivatives  and  moments  in  the  budgets  can  be  computed 
with  confidence.  This  is  particularly  important  as  the  current  computations  employ  novel 
LBM  approach  rather  than  the  tried  and  tested  Navicr-Stokcs  schemes.  Figure  34  shows  the 
comparison  of  time  derivative  of  Reynolds-strcss  (U1U2)  computed  from  two  independent  means. 
One  computation  comes  from  summing  the  terms  on  the  right  hand  side  of  Reynolds  stress 
equation.  The  second  estimate  of  the  time-derivative  comes  from  the  finite  difference  operation 
on  Reynolds-strcss: 

d  ( UiUj )  ( UiUj(t  +  dt))  —  ( UiUj(t  —  dt)) 

dt  “  2 Jt  * 

The  results  from  these  two  different  computations  are  indistinguishable.  Further  tests  were 
performed  to  ensure  the  fidelity  of  other  budget-term  calculations  (figures  not  shown).  The 
results  clearly  indicate  that  all  terms  in  the  Reynolds-strcss  budget  equation  arc  adequately 
accurate. 

Figure  35  shows  different  terms  in  the  budget  of  ( u\u2 )  evolution.  The  magnitudes  of  the 
various  terms  diminish  with  time  as  k  decays  from  its  initial  value.  As  is  to  be  expected 
in  the  case  of  off-diagonal  stress,  dissipation  is  quite  small  at  all  times.  The  pressure-strain 
distribution  always  counters  production  and  the  difference  between  the  two  drives  the  evolution. 
This  balance  is  very  similar  to  that  seen  in  the  constant  shear  case  (Rogers  1986),  with  the 
exception  that  production  and  pressure-strain  distribution  arc  now  cyclic. 

The  budget  of  (u\U\)  evolution  is  shown  in  Fig.  36.  As  in  the  case  of  constant  shear, 
dissipation  dominates  the  very  early  evolution  (5max£  <  1)  as  the  Reynolds  stress  value  falls 
rapidly.  Although  zero  initially,  production  becomes  the  most  dominant  term  soon  (5max<  >1). 
Interestingly,  the  peak  value  of  normalized  production  in  the  time-varying  and  constant  shear 
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cases  occur  at  the  same  normalized  time:  Smaxt  =  2.  Pressure-strain  redistribution  almost 
always  counters  the  effect  of  production,  but  there  are  rare  instances  when  both  the  processes  are 
of  the  same  sign.  It  remains  to  be  seen  if  the  current  pressure-strain  redistribution  closures  can 
capture  this  behavior  as  most  models  assume  that  redistribution  always  counteracts  production. 
At  the  intermediate  stages  (1  <  Smax£  <  10),  the  difference  between  production  and  the  sum  of 
dissipation  and  pressure-strain  redistribution  drives  the  evolution.  At  later  stages,  dissipation 
becomes  negligible  (in  this  decaying  ease). 

The  time  dependence  of  budgets  of  {U2U2)  and  (u3u3)  arc  quite  similar  and  latter  budget  is 
shown  in  Fig.  37.  For  both  components,  production  is  absent.  Dissipation  dominates  the  initial 
time  ( Smaxt  <  3)  development  but  quickly  goes  to  zero  thereafter.  The  later  development  is 
mainly  due  to  pressure-strain  redistribution.  The  production  of  ( U\U\ )  is  the  only  source  of 
kinetic  energy  in  this  flow.  We  see  that  the  magnitudes  of  production,  the  dissipation  and 
redistribution  decrease  with  time.  The  increase  of  the  magnitude  of  P/e  in  Fig.  28  is  caused 
by  the  faster  decrease  of  e,  rather  than  the  increase  of  P. 

Overall,  it  is  clear  that  pressure  redistribution  plays  a  key  role  in  determining  the  evolution 
of  various  normal  components.  Calculations  were  performed  to  confirm  that  the  trace  of  the 
pressure-strain  redistribution  was  indeed  close  to  zero.  Analysis  of  the  three  components  reveals 
that  Ru  and  R22  arc  nearly  in  phase,  while  the  phase  of  i?33  is  shifted  nearly  by  7r.  This  means 
that  energy  is  removed  from  (u\U\)  and  (1*21*2)  and  added  to  (u3u3).  This  is  different  from 
what  is  observed  in  constant  shear  homogeneous  flow.  In  the  constant  shear  case,  pressure- 
strain  redistribution  removes  energy  from  (uiUi)  and  augments  (112U2)  and  (u3u3).  Closer 
examination  of  the  constant  shear  budgets  of  Rogers  (1986)  reveals  that  the  rapid  pressure- 
strain  redistribution  docs  remove  energy  from  both  ( U\U\ )  and  (1*2^2)  and  augments  (u3u3).  It 
is  the  slow  pressure-strain  redistribution  that  is  responsible  for  the  augmentation  of  the  (1*21*2) 
component.  In  the  present  time-varying  shear  ease,  we  expect  the  non-linear  slow  term  to  be 
less  dominant  at  later  times  as  the  turbulence  decays  quickly. 

3.1.6  Limiting  high-frequency  behavior 

As  mentioned  in  the  Introduction,  the  present  flow  shares  some  similarities  with  homogeneous 
shear  flow  in  a  rotating  reference  frame.  For  example,  in  both  eases  turbulence  grows  when  the 
frequency  (of  rotation  or  periodic  variation)  is  small  and  decays  when  the  frequency  is  large. 
However,  there  are  very  significant  differences  in  the  turbulence  dynamics  in  the  two  cases.  In 
the  case  of  homogeneous  shear  flow  in  a  rotating  reference  frame,  turbulence  tends  to  a  two- 
componential  limit  at  large  rotation  rates  in  keeping  with  the  Taylor-Proudman  theorem.  Two- 
dimensionalization  inhibits  cascade  (no  vortex  stretching)  and  diminishes  production  leading  to 
viscous  decay  of  the  velocity  field  at  all  scales.  In  the  present  case,  on  the  other  hand,  turbulence 
is  isotropic  at  large  frequencies  of  shear  variation.  The  shear  varies  too  rapidly  to  induce 
anisotropy  and  produce  kinetic  energy.  We  anticipate  that  at  high  frequencies,  turbulence 
subject  to  rapidly  varying  shear  will  be  more  like  decaying  isotropic  turbulence  as  production 
will  not  be  significant.  To  establish  that  this  is  indeed  the  case,  we  will  now  compare  u;/Smax=10 
case  with  decaying  isotropic  turbulence. 

In  the  decaying  isotropic  turbulence  simulation,  all  the  initial  conditions  are  the  same  as  for 
homogeneous  shear  flow,  except  that  wall  boundary  condition  is  replaced  by  periodic  boundary 
condition  and  there  is  no  mean  shear  in  the  flow.  Figure  38  shows  the  evolution  of  k  and  e 
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for  homogeneous  shear  and  isotropic  decaying  cases.  For  consistency,  time  is  normalized  with 
Smnx  >  although  there  is  no  shear  in  isotropic  decaying  turbulence.  In  homogeneous  shear  ease, 
at  initial  stage,  k  decays  a  little  bit  slower  than  in  isotropic  decaying  case.  But  over  all,  the 
decay  rates  arc  almost  same  for  both  eases.  A  similar  examination  of  the  dissipation  rate 
demonstrates  that  the  evolution  in  the  two  eases  arc  almost  indistinguishable.  The  standard 
power-law  decay  is  recovered  in  both  cases.  It  must  be  pointed  out  that  there  are  about  15 
shear  cycles  in  the  period  over  which  the  comparison  is  made. 

Next  we  compare  the  evolution  of  anisotropy  in  the  two  eases  in  Fig.  39.  The  small  level 
of  anisotropy  seen  in  the  initial  condition  is  due  to  finite  statistical  sample  size.  The  time- 
development  of  the  diagonal  anisotropies  in  the  decaying  and  high-frequency  forcing  eases  arc 
quite  similar,  with  the  latter  exhibiting  small-amplitude  oscillations  about  the  former.  With 
regard  to  612,  the  difference  between  the  two  cases  is  quite  discernible.  The  negative  bias  of 
b\2  in  the  forced  case  has  already  been  explained  above.  The  most  significant  features  are:  the 
cycle-averaged  612  in  the  forced  case  is  quite  different  from  the  &i2  in  the  isotropic  decaying  ease; 
and  the  amplitude  of  oscillation  is  quite  sizable.  Despite  the  large  amplitude  of  612  variation, 
the  production  is  very  small  due  to  the  fact  that  stress  and  strain  are  out  of  phase  by  7r/2. 
Thus,  the  diagonal  anisotropies,  kinetic  energy  and  dissipation  are  not  very  different  from  those 
in  the  unforced  decaying  ease. 

Thus,  the  limiting  turbulence  behavior  in  homogeneous  shear  flow  in  a  rotating  reference 
frame  is  dictated  by  Taylor-Proudman  theorem  and  characterized  by  two-dimensional  fluctua¬ 
tions.  In  the  periodic  shear  ease,  the  limiting  behavior  is  characterized  by  statistically  isotropic 
fluctuations  and  turbulence  decays  as  in  unforced  case. 

Summary  of  Study  2.  The  RDT  results  reveal  that  the  turbulence  response  to  this  type 
of  rotating  shear  is  somewhat  unexpected.  We  observe  the  following  responses  as  a  function  of 
forcing  frequency:  (i)  At  low  forcing  frequencies  (0  <  A /S  <  0.14),  initially  isotropic  turbu¬ 
lence  rapidly  evolves  to  one-component  state  and  kinetic  energy  asymptotes  to  a  value  smaller 
than  the  initial  level,  (ii)  At  intermediate  frequencies  (0.14  <  X/S  <  0.20),  the  turbulence 
kinetic  energy  exhibits  sporadic  spurts  of  growth,  (iii)  At  larger  frequencies,  production  and 
pressure-strain  distribution  perturb  the  field  only  slightly  from  its  isotropic  initial  state.  The 
sporadic  growth  observed  at  the  intermediate  forcing  frequencies  is  quite  unexpected.  A  quali¬ 
tative  explanation  for  the  observed  behavior  is  presented  in  the  referenced  paper.  It  is  argued 
that  the  rapid  pressure-strain  correlation  is  responsible  for  this  somewhat  anomalous  behavior. 
Although  the  physics  of  the  observed  behavior  is  quite  complicated,  we  demonstrate  that  by 
simply  accounting  for  the  phase  lag  between  stress  and  strain,  the  observed  behavior  can  be 
approximately  reproduced.  A  more  quantitative  explanation  will  require  precise  closure  model 
for  the  rapid  pressure-strain  correlation.  None  of  the  current  models  come  close  to  predicting 
the  behavior  observed  in  this  paper.  Development  of  closure  models  that  can  quantitatively 
capture  the  observed  physics  is  the  next  step. 
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Figure  16:  Computational  domain.  The  statistics  arc  obtained  in  the  core  region  of  computa¬ 
tional  domain  which  is  128  x  64  x  128. 
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Figure  17:  The  profile  of  mean  velocity  at  St  =  12  and  the  distribution  of  k(y )  at  initial 
stage  and  St  —  12.  The  wall  effects  are  limited  to  regions  near  the  walls.  In  the  inner  region, 
mean  velocity  agree  very  well  with  the  ideal  mean  velocity  and  homogeneity  of  turbulence  is 
maintained.  , 
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Figure  18:  Evolution  of  normal  Reynolds  stress  in  homogeneous  turbulent,  shear  flow  (constant 
shear).  NS-DNS  data  are  from  Rogers  (1986)  with  initial  values  of  Rex  ~  35  and  5*  =  1.2. 
Experimental  data  are  from  Tavoularis  and  Corrsin  (1981a)  with  ReA  =  284  at  St  =  8.6. 


Figure  19:  The  results  of  612  obtained  from  DNS  using  LBM  and  Navier-Stokcs  equations  by 
Jacobitz,  Sarkar  &  Van  Atta  (1997)  with  the  initial  values  of  RcA  =  44.72  and  Sk/e  —  2.0. 
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Figure  20:  The  profiles  of  U(y)  and  k(y)  at  Smaxt  =  34.55,  37.70,  and  40.81  for  ui/Smax- 0.5. 
Ideal  vs.  computed. 
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Figure  21:  Planar- averaged  Reynolds  stresses  along  y-dircction  at  Smaxt  =  34.55. 
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Figure  22:  Planar-avcragcd  Pressure-strain  correlations  along  y-dircction  at  Smaxt  =  34.55. 
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Figure  23:  Planar-avcragcd  components  of  dissipation  tensor  along  y-dircction  at  Smaxf  =  34.55. 


Figure  24:  Evolution  of  k  in  various  u/ Smnx  cases. 
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Figure  25:  Evolution  of  k  and  e  in  u/Smax  =  0.5  ease. 


Figure  26:  Evolution  of  k  and  e  in  u)/Smax  =  1.0  ease. 
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S/S 


Figure  27:  Evolution  of  production  over  dissipation  ratio  (P/e)  in  u/Smax  =  0.5  case. 


Figure  28:  Evolution  of  production  over  dissipation  ratio  (P/e)  in  uj/Smax  =  1.0  case. 
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Figure  29:  Evolution  of  b\2  in  u/Smax  =  0.5  case. 


Figure  30:  Evolution  of  bi2  in  u)/Smax  =  1.0  ease. 
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Figure  31:  Evolution  of  6n  and  b22  in  uj/Smax=  0.5  case. 


34 


S/S 


b„ 


Figure  32:  Evolution  of  6n  and  b2 2  in  u)/Smax=  0.65  case. 


Figure  33:  Evolution  of  bu  and  b22  in  u>/Smax—  10  case. 
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Figure  34:  Comparison  of  d  (iiiUi)  /dt  computed  from  two  different  methods  in  u/Smax  =  10 
case. 


Figure  35:  Temporal  variation  of  (uiu2)  budget  in  uj/Smax—  10  case. 
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Figure  36:  Early-time  variation  of  (ujU])  budget  in  u>/Smax  =  1.0  ease. 


Figure  37:  Early-time  variation  of  (u3u3)  budget  at  earlier  stage  in  u >/Smax  =  1.0  ease. 
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Figure  38:  The  evolutions  of  k  and  e  in  isotropic  decaying  turbulence  and  homogeneous  shear 
turbulence  with  u/Smax  =  10. 


Figure  39:  Evolution  of  anisotropy  tensor  in  isotropic  decaying  turbulence  and  homogeneous 
shear  turbulence  with  u;/Smax  =  10. 
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4  Unfinished  problems. 

We  now  describe  the  unfinished  problems  that  were  initiated  during  the  last,  year  of  the  work. 

4.1  Near  wall  unsteadiness  issues 

Fan  &  Lakshminarayana  [1]  and  Mankbadi  &  Mobark  [2]  consider  unsteady  boundary  layers 
where  either  the  wall  or  the  free  stream  fluctuates.  The  flow  is  characterized  by  a  near-wall 
flow  that  fluctuates  in  phase  with  the  wall  and  an  outer  region  with  a  substantial  phase  shift. 

An  instantaneous  log-law  docs  not  in  general  exist  and  near-wall  formulations  based  on 
the  log-law  and  the  equilibrium  assumption  arc  not  appropriate  according  to  Fan  &:  Laksh¬ 
minarayana.  That  disqualifies  algebraic  eddy-viscosity  models  and  two-equation  models  with 
log-law  boundary  conditions  or  with  near  wall  damping  functions  based  on  the  wall  friction 
(such  as  y+).  Moreover  the  existing  low  -Re  K  models  are  developed  for  steady  flows  and  do 
not  give  satisfactory  results  for  unsteady  boundary  layers. 

Requirements  on  RANS  models  for  highly  unsteady  turbulence.  When  the  frequency 
of  large-scale  unsteadiness  increases  and  approaches  the  turbulence  characteristic  frequency, 
turbulence  is  profoundly  affected.  That  means  that  non-equilibrium  effects  become  more  pro¬ 
nounced.  The  very  extreme  kind  of  unsteady  turbulence  is  rapidly  distorted  turbulence  which 
has  been  thoroughly  investigated  in  the  first  part  of  the  study.  It  is  reasonable  to  expect  that 
turbulence  models  that  arc  suited  for  non-equilibrium  turbulence  will  perform  well  in  unsteady 
turbulent  flows  as  well.  Some  discussion  in  this  regard  can  be  found  in  Mankbadi  &:  Mobark 
[2].  In  keeping  with  the  objectives  of  current  research,  we  identify  the  following  physical  flow 
features  as  being  crucial  for  successful  RANS  computation  of  unsteady  wall-bounded  flows: 

1.  The  unsteady  wall  shear  stress  is  out  of  phase  with  the  flow  away  from  the  wall.  Explicit 
corrections  of  the  near  wall  influence  cannot  be  based  on  the  wall  friction  velocity  uT  . 
The  local  turbulent  velocity  scale  \[K  may  be  an  alternative  and  near-wall  functions  may 
be  based  on  Ry  =  \J~Kylv  or  Rt  =  K2 /vc.. 

2.  The  near- wall  unsteady  region  decreases  in  size  for  increasing  frequencies  and  could  be  of 
the  order  of  the  viscous  sub-layer.  Correct  and  consistent  near-wall  asymptotic  behavior 
is  thus  important  for  correct  prediction  of  the  mean  wall  shear  stress.  It  is  also  important 
that  the  turbulence  profiles  are  well  reproduced  through  the  buffer  layer  and  thus  models 
such  as  the  Jones  &  Launder  and  the  Wilcox  standard  K  —  u  models  are  not  well  suited 
due  to  the  quite  bad  reproduction  of  the  peak  in  the  K-profile. 

3.  The  prediction  of  boundary  layer  separation  is  important  for  capturing  the  unsteady  flow 
around  transonic  airfoils.  The  eddy-viscosity  assumption  fails  to  correctly  predict  the 
production  of  the  turbulent  kinetic  energy.  That  can  instead  be  obtained  by  some  ad  hoc 
limitation  of  the  production  or  eddy-viscosity  but  the  same  effect  is  obtained  by  using 
models  based  on  full  RST  models. 

4.  When  the  frequency  of  the  unsteadiness  increases,  the  capability  of  the  model  to  predict 
nonequilibrium  turbulence  increases  in  importance.  That  is  not  very  well  fulfilled  by 
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standard  cddy-viscosity  two-cquation  models,  but  somewhat,  improved  in  c.g.  the  Menter 
SST  K  -  u  model.  Also  here  models  based  on  full  RST  models  arc  better  suited  for 
describing  such  flows. 

Objective  and  approach.  The  main  objective  of  this  research  is  to  seek  physics-based  mod¬ 
ifications  to  current  two-equation  and  second  moment  closure  (SMC)  models  to  improve  their 
performance  in  unsteady  flows.  In  such  flows,  the  quasi-steady  assumption  will  be  inadequate. 
In  this  research  period,  we  explored  the  effects  of  two  important  changes  to  the  standard  tur¬ 
bulence  model  paradigm.  First  and  foremost,  we  recognize  that  spectral  transfer  rate  and 
dissipation  are  not  equal  to  one  another  in  statistically  unsteady  turbulent  flows.  Given  the 
nature  of  spectral  energy  cascade,  it  can  be  easily  argued  that  dissipation  lags  behind  cascade 
rate  when  the  energy-containing  range  of  the  spectrum  undergoes  rapid  changes.  This  has 
significant  impact  on  the  turbulence  length  scale  equation  and  we  propose  a  2-equation  length- 
scale  model  to  account  for  these  effects.  This  will  have  an  impact  on  the  current  two-equation 
(which  will  become  3-cquation)  and  SMC  closures.  An  important  inference  from  the  above 
list  of  desirable  features  mentioned  above  and  our  previous  DNS  (direct  numerical  simulations) 
and  RDT  (rapid  distortion  theory)  studies  is  that  there  could  be  a  significant  phase  difference 
between  the  imposed  mean  velocity  gradient  (mean  strain)  and  Reynolds  stress.  This  is  due  to 
the  fact  that  the  Reynolds  stress  depends  on  the  strain-rate  history  rather  than  the  local  mean 
strain  rate.  Thus,  in  our  second  investigation,  we  study  the  effect  of  introducing  a  phase  lag 
between  stress  and  strain  in  the  Boussinessq  constitutive  relation.  This  work  will  only  affect 
current  two-cquation  models.  This  phase-lag  can  appear  intrinsically  in  SMC  models. 

4.2  Two-equation  length-scale  model 

The  standard  K  —  e  model  is  written  as 

P-e  +  DK, 

c2 

P> :  ~  i  (19) 

where  Pt  is  the  production  of  the  dissipation  e.  In  the  standard  K  -  e  approach,  the  turbu¬ 
lence  spectra  is  assumed  to  be  in  equilibrium,  and,  thus,  P(  may  be  directly  related  to  P,  the 
production  of  K 

eP 

P*  =  ~K  (2°) 

However,  in  non-equilibrium  flows  the  e  equation  should  not  immediately  respond  to  a 
sudden  change  in  production  of  K.  Therefore,  we  propose  to  add  relaxation  to  the  production 
of  e  by  adding  an  additional  transport  equation  for  s  =  Pf. 

The  resulting  three-equation  model  may  be  written  as 
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Dt 

40 


=  P  -  f  +  Dk  , 

e2 

=  s  —  Ct2  —  +  D( , 


DK 

~Dt 

De 

~Dt 


Ds 

Dt 


c 


(21) 


K 


cP 

1~K~ 


+  Ds . 


At  this  stage,  let  us  neglect  the  diffusion  of  s  (Ds  -  0).  The  relaxation  time  scale  is 
proportional  to  the  turbulent  time  scale  K/(.  with  the  relaxation  proportionality  constant  Cs , 
being  of  order  unity.  This  model  will  have  the  correct  asymptotic  behavior  in  homogeneous 
decaying  turbulence.  However,  the  asymptotic  state  in  homogeneous  developing  flow  will  not 
be  obtained.  The  reason  is  that  the  production  of  e  will  not  be  constant  at  equilibrium.  In  order 
to  obtain  the  correct  asymptotic  state,  we  may  let  s  represent  the  production  of  r.  normalized 
by  K ,  and,  thus,  s  becomes  constant  at  the  asymptotic  state.  The  resulting  three-equation 
model  then  becomes 
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(22) 


In  order  to  solve  these  equations  in  wall  bounded  flows  we  need  some  near-wall  treatment.  The 
Wilcox  standard  K  —  u  model  is  one  of  the  preferred  models,  since  it  does  not  rely  on  friction 
velocity  scaling.  Let  us  transform  equations  (22)  to  a  K  -  u>  -  s  form,  where  u  is  defined  as 


u  = 


e 

Wk' 


(23) 


The  transformed  model  reads 
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The  Wilcox  K  -  u  model  coefficients  are  related  to  the  K  -  e  coefficients  by  the  following: 


0'  =  =  0.09  ^  =  0.09 , 

a  =  (Cel  -  1)  =  0.556  ^  =  1.556 , 

0  =  0'(Cl2~  1)  =  0.075  ^  Cl2  =  1.883  .  (25) 

The  author  would  like  to  further  pursue  this  work  in  the  future.  (The  three-equation  model 
was  developed  in  collaboration  with  Dr.  Stefan  Wallin  of  FOI/KTH,  Sweden.) 


4.3  Phase-lagged  Boussinessq  Approximation 

As  demonstrated  in  our  previous  papers  (published  previously  under  this  project),  in  highly 
unsteady  flows,  Reynolds  stress  lags  behind  mean-flow  strain.  This  is  due  to  the  fact  that  the 
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stress  depends  on  the  strain-rate  history  rendering  the  standard  Boussinessq  approximation 
invalid.  We  are  currently  investigating  a  model  described  by  the  following  equations: 


dk 

dt 

de 

dt 

(u,Uj) 


p 


P-e, 

P  e  e2 

Ccl—  -ce2-, 

k  k 

lk8x]-2C^Sl] 

2  k 2 

—  k  Sjj  —  2Cll—Sij(t  —  p) , 

=  2 i'TSij(t)Sij{t  -  p) . 


(26) 


In  future  work,  the  author  will  demonstrate  this  model  possesses  many  important  features  that 
have  not  yet  been  fully  exploited  in  current  literature. 
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